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Quantum degenerate gases trapped in optical lattices are ideal testbeds for fundamental physics because 
these systems are tunable, well characterized, and isolated from the environment. Controlled disorder can be 
introduced to explore suppression of quantum diffusion in the absence of conventional dephasing mechanisms 
such as phonons, which are unavoidable in experiments on electronic solids. Recent experiments use transport 
of degenerate Fermi gases in optical lattices (Kondov et al. Phys. Rev. Lett. 114, 083002 (2015)) to probe a 
particularly extreme regime of strong interaction in what can be modeled as an Anderson-Hubbard model. These 
experiments find evidence for an intriguing insulating phase where quantum diffusion is completely suppressed 
by strong disorder. Quantitative interpretation of these experiments remains an open problem that requires 
inclusion of non-zero entropy, strong interaction, and trapping. We argue that the suppression of transport can be 
thought of as localization of Hubbard-band quasiparticles. We construct a theory of transport of Hubbard-band 
quasiparticles tailored to trapped optical lattice experiments. We compare the theory directly with center-of- 
mass transport experiments of Kondov et al. with no fitting parameters. The close agreement between theory 
and experiments shows that the suppression of transport is only partly due to finite entropy effects. We argue that 
the complete suppression of transport is consistent with Anderson localization of Hubbard-band quasiparticles. 
The combination of our theoretical framework and optical lattice experiments offers an important platform for 
studying localization in isolated many-body quantum systems. 

PACS numbers: 03.75.Ss, 67.85.-d 


I. INTRODUCTION 


Understanding the motion of a quantum particle in an oth¬ 
erwise isolated lattice under the influence of an applied field 
is central to our understanding of conductivity in electronic 
solids. The theory of Anderson localization [1,2] predicts 
that quantum diffusion of a single particle can fail in a dis¬ 
ordered lattice. Above a critical disorder strength, for which 
the mobility edge encompasses all states participating in trans¬ 
port [3, 4], strong interference forbids quantum diffusion. An¬ 
derson’s mechanism of localization was first discussed in the 
context of as a simplified model designed to treat the propa¬ 
gation of highly excited states of nuclear spin systems but is 
has since been applied to a wide variety of other systems [2], 
including quantum degenerate atomic gases [5-10], Disorder- 
induced localization is also believed to play a key role in 
metal-insulator transitions in a wide-range of materials [2-4], 

Subsequent theoretical studies of Anderson localization 
found that inclusion of realistic effects, specifically inter¬ 
particle interactions and non-zero temperature [2, 11- 5], 
pose prominent problems. The competition between Ander¬ 
son localization and strong interaction effects have been stud¬ 
ied with a variety of methods, e.g., quantum Monte Carlo [ 1 6], 
dynamical mean field theory [14, 17, 18], and related quantum 
cluster methods [19], Refs. [14] and [18], for example, found 
a correlated Anderson insulator ground state for large disor¬ 
der strengths indicating that Anderson localization persists in 
a strongly interacting limit. A more complete understanding 
of the interplay of strong inter-particle interactions and disor¬ 
der is urgently needed to enhance our knowledge of strongly 
correlated materials such as high-temperature superconduc¬ 
tors. 

Related work by Basko et al. [20] has triggered consider¬ 



FIG. 1. (Color online) Schematic showing disordered lattice sites in 
a parabolic trapping potential. The site coloring represents a dense 
core that gives way to zero density at the edges. The system stud¬ 
ied here can be thought of as a strongly interacting high temperature 
paramagnet with a density less than one at the center. An applied 
shift of the external trapping potential along the ^-direction for a time 
r = Tp forces center-of-mass motion along the x direction only if 
the atoms are mobile. Tp is chosen to be short on the time scale of 
the inverse trapping frequency. 


able interest in the interplay between interactions, tempera¬ 
ture, and Anderson localization. Their work indicates that a 
correlated Anderson insulator is stable at non-zero tempera¬ 
tures and corresponds to a many-body localized state. This is 
surprising because one might expect that interactions lead to 
dephasing effects that mimic the effects of heat and particle 
number reservoirs [21] that are known to lead to conduction 
via variable range hopping in certain solids [4]. Interactions 
would be expected to lead to effective reservoirs even in the 
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absence of an explicit reservoir. But Ref. [20] argues, surpris¬ 
ingly, that interactions allow a correlated Anderson insulator 
to survive up to a characteristic temperature. 

Quantum degenerate gases of atoms trapped in optical lat¬ 
tices offer a controlled arena to study the interplay of inter¬ 
actions, disorder, and thermal effects [22-25] because they 
are, to an excellent approximation, isolated. The entropy 
per particle, controlled via cooling in a parabolic trap, deter¬ 
mines an equilibrium temperature when the lattice is turned 
on, since atomic gases thermalize through inter-particle inter¬ 
actions [26, 27], As a result of their isolation, quantum de¬ 
generate Fermi gases in optical lattices exhibit quantum diffu¬ 
sion (see, e.g.. Ref [28]), even though their temperatures are 
a significant fraction of the Fermi energy. This offers a use¬ 
ful regime to study because isolated systems can, in princi¬ 
ple, exhibit many-body localization even at high temperature 
[29]. Furthermore, optical lattice experiments are well char¬ 
acterized [24, 30]: interaction strength, lattice depth, entropy, 
density, and other parameters are all known and tunable. The 
impact of disorder can therefore be studied independently of 
conventional dephasing phenomena arising from contact with 
reservoirs [31-38], 

Recent experimental work [39] has investigated interacting 
fermions confined in a cubic optical lattice to study the influ¬ 
ence of quenched speckle disorder on center-of-mass trans¬ 
port, Fig. 1. This system is accurately described using the 
Anderson-Hubbard model. They find intriguing insulating be¬ 
havior above a critical disorder strength that agrees qualita¬ 
tively with the predictions for many-body localization in the 
weakly interacting regime [20]. Control over the lattice poten¬ 
tial depth and disorder strength was shown to lead to a regime 
where known types of insulating behavior can be excluded. 
For example, it was demonstrated that the insulating regime 
occurs for disorder strengths well below the classical percola¬ 
tion threshold [39]. Furthermore, the system was made dilute 
enough to avoid forming Mott [40-42] and band insulators. 
The regime they explored can be thought of as a strongly inter¬ 
acting Hubbard paramagnet with a temperature well below the 
bandwidth. The insulating behavior of the isolated-strongly 
interacting particles in these experiments [39] is therefore a 
highly non-trivial probe of localization in many-body quan¬ 
tum states. 

In this article, we provide, to our knowledge, a new per¬ 
spective on disorder-induced localization in the Anderson- 
Hubbard model and the measurements in Ref. [39]. We es¬ 
tablish a connection to Hubbard-band quasiparticles using a 
direct comparison between theory and experiment with no fit¬ 
ting parameters. This approach enables us to treat the strongly 
interacting limit, which was not possible using the perturba¬ 
tion theory employed in Ref. [20], We derive the equations 
of motion for the Anderson-Hubbard model in the paramag¬ 
netic regime while taking into account all important exper¬ 
imental aspects, particularly trapping and finite entropy ef¬ 
fects. We show that the equations of motion derived here 
to include a trap reduce to the Hubbard-I [43] approximation 
normally considered in the translationally invariant limit. This 
demonstrates that our equations of motion quantitatively cap¬ 
ture the dynamics of Hubbard-band quasiparticles in a trap. 


The Hubbard-I approximation is a strong coupling approxi¬ 
mation that becomes exact in the limit of strong interactions 
and high entropy. Our formalism is therefore directly applica¬ 
ble to the experiments of Ref. [39], 

We use parameters taken from Ref. [39] to effectively repli¬ 
cate the experiment numerically and find evidence for a quasi¬ 
particle mobility edge. We find that at low disorder strengths 
the quasiparticles propagate in the lattice under an applied 
force, i.e., they have non-zero mobility. We also study the 
result of increasing disorder. At large disorder strengths we 
identify a transition to an insulator through the absence of 
center-of-mass motion. A direct comparison between the¬ 
ory and experiment shows good agreement. We argue that 
the insulating behavior observed in Ref. [39] is consistent 
with Anderson localization of Hubbard-band quasiparticles. 
To our knowledge, disorder-induced localization in the Hub¬ 
bard model has not been previously understood using this 
approach, which is complementary to other methods, e.g., 
perturbative theory [20] and dynamical mean field theory 
[14,17,18], 

We begin in Sec. II by defining the model used to simu¬ 
late the experiments of Ref. [39] and all necessary parame¬ 
ters. Here we also define the center-of-mass velocity as the 
key observable. In Sec. Ill we then derive the equations of 
motion in the paramagnetic regime. Sec. IV then shows that 
the equations of motion reduce to the Hubbard-I approxima¬ 
tion [43] that was originally used to define Hubbard-band 
quasiparticles. Here we also show that, in a strong cou¬ 
pling limit, Hubbard-band quasiparticles obey an effective 
Anderson model of non-interacting quasiparticles. Sec. V 
then defines the approximations used in constructing the ini¬ 
tial state that is propagated using the equations of motion. 
Sec. VI points out an important feature of the initial states 
used in these experiments. We find that increasing disorder 
at fixed entropy effectively raises the system temperature to 
at most B/ 3, where B is the bandwidth. Even though this 
heating keeps the temperature well below the bandwidth, it 
is nonetheless an important aspect of these experiments that 
must be included to make a quantitative comparison with the¬ 
ory. Sec. VII presents our central results. Here we directly 
compare numerical solutions of the equations of motion with 
experiments. We find that low disorder allows the Hubbard- 
band quasiparticles to propagate in the trap. But we find a crit¬ 
ical disorder strength above which all transport is suppressed. 
We conclude in Sec. VIII by interpreting the results presented 
here as evidence for the Anderson localization of Hubbard- 
band quasiparticles. 


H. MODEL AND PARAMETER REGIMES 

We study the dynamics of an equal population of two- 
component fermions in a cubic optical lattice in the presence 
of spatial disorder. For deep lattices we assume that all N par¬ 
ticles reside in the lowest Bloch band. In this limit the single¬ 
band Anderson-Hubbard model is an excellent approximation 
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[24, 39]: 

h ah = T j,i' c ],v c j',o + U Y1 n rt n U + Vo n i (!) 

j,j',cr j j 


Here cj CT creates a fermion of spin a £f, ]. at a site Rj, U is a 
repulsive interaction, and rij = nj y ^+rij^ is the number oper¬ 
ator. The matrix elements Tjji = ir) R R ,4 ,5 are written in 
terms of a Kronecker delta, 5, that enforces a hopping energy 
t between nearest neighbor sites (S is a nearest neighbor bond 
vector.) 

The last term in Eq. (1) includes spatially inhomogeneous 
perturbations to the chemical potential. We define: 

muj 2 a 2 R 2 

Hj = —/j ,0 H- 2 ^ e i P VpR-j 1 x, (2) 

where fj, 0 is the average chemical potential, m is the atomic 
mass, oj is the trapping frequency that parameterizes the ex¬ 
ternal confinement, a is the lattice spacing determined by the 
optical lattice laser wavelength, ej denotes spatially random 
disorder, and Vp is a pulse strength that is switched on for a 
time t = Tp to effectively shift the trap center. 

Vp acts as the analogue of a weak electric field used to 
drive transport along the x-direction, see Fig. 1. At long 
times a single particle with no disorder will oscillate in the 
trap. But we consider pulse times that are short with respect 
to the inverse trapping frequency to focus on the linear re¬ 
sponse regime probed in Ref. [39], At these short times, the 
center-of-mass velocity is unidirectional and, in the absence 
of disorder, increases linearly with Vp. 

We consider two distinct distributions of site disorder. In 
the experiments of Ref. [39] the speckle potential used to es¬ 
tablish a disordered optical lattice creates an exponential prob¬ 
ability distribution function for the on-site energies [30]: 


P E (e) 


g— e/As 

Ae 


(3) 


where Ap is the strength of the exponentially distributed dis¬ 
order assuming e > 0 (this is accurate to within 10 % of the 
disorder strength used in Ref. [39]). We also consider a uni¬ 
form (boxed) disorder probability distribution function for the 
on-site energies e :l : 


Pu(e) 


e(A{//2 — |e|) 

A u 


(4) 


where 0 is the Heaviside step function and Ap parameter¬ 
izes the strength of the uniformly distributed disorder. Pp in¬ 
troduces behavior that is distinct from more common models 
with Pu because changing Ap at fixed N forces // (l to change. 
This is in contrast to changes in Ap which leaves fio constant 
at fixed N. 

Eq. (1) quantitatively captures the essential properties of 
the experiments in Ref. [39], We ignore disorder in t and U 
that was shown [30] to be narrowly Lorentzian distributed. In 
what follows, we find that we are able to make quantitative 
comparison with experiment even while excluding the disor¬ 
der in t and U. We will return to this point in Sec. IV. 


TABLE I. Table of parameters used in the experiments of Ref. [39]. 
Here the recoil energy is Er = h 2 / ( 8 ma 2 ) and the atomic species 
is 40 K. 


Lattice Depth 

V L 

6 Er 

7Er 

Trap Frequency 

UJ 

110 x 2ttHz 

114 x 2ttHz 

Lattice Spacing 

a 

391.1 nm 

391.1 nm 

Number of Particles 

N 

47100 ± 6500 

48700 ±1900 

Entropy per Particle 

S/N 

1.51 ± 0.18fc B 

1.6 ±0.17 k B 

Hopping 

t 

0.0509 Er 

0.0395 Er 

Interaction 

U 

0.304 Er 

0.355 Er 

Relative Strength 

U/t 

5.97 

8.98 

Disorder Strength 

A E 

0-2 Er 

0-2 Er 

Pulse Time 

rp 

2 ms 

2 ms 

Pulse Strength 

Vp 

0.011 Er 

0.011 Er 


The experiments proceed by trapping a fixed number of par¬ 
ticles with a fixed entropy, S. The entropy and all other nec¬ 
essary model parameters were determined in Ref. [39] and are 
shown in Table I. We focus on the two lattice depths with high 
U, where U/t ss 6 and U/t ss 9 for 6 Ep and 7 Ep, respec¬ 
tively, which allows for a strong coupling approximation that 
becomes exact in the limit U/t — » 00. Table I leaves no fit¬ 
ting parameters in using approximate solutions of Eq. (1) to 
compare with the experiments of Ref. [39]. 

We will show that the entropies reported in Table I im¬ 
ply temperatures that are well above the Neel temperature, 
~ t 2 /U [42, 44-48], The experimentally relevant tempera¬ 
ture regimes are above the hopping but below the bandwidth. 
Our central set of approximations in studying Eq. (1) can be 
summarized by: 

t < U 

t <k B T< 12 1 (5) 

where the first inequality assumes that we focus on the high 
lattice depth data of Ref. [39], and the second inequality im¬ 
plies that high temperature limits are valid approximations. 
Sec. V shows that the initial state for the parameters defined 
by Table I can be thought of as a dilute ((n) < 1) high tem¬ 
perature paramagnet. We will therefore focus our study to 
strongly interacting paramagnetic regimes. 

To make contact with experimental results presented in 
Ref. [39] we study the dynamics of the center of mass. The 
time dependent center-of-mass velocity in particular: 

V c .o.M.(r) = ^R i «n J -)) D ( 6 ) 

0 

was inferred from time of flight images [39], Here ((■■■))d 
indicates disorder averaging of expectation values and r de¬ 
notes time. In the following we find that disorder averaging 
over 25-50 realizations is sufficient to reach convergence in 
our numerical simulations. We will use Eq. ( 6 ) to compute 
the center-of-mass velocity along the direction of the applied 
pulse after a time r = rp: 

Vc.o.m. = x ■ Vc.o.m.(tp)- (7) 
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This quantity is akin to measures of mobility in solids. For 
example, in the Drude model of electrical conductivity, Vc.o.m. 
is proportional to the electron mobility when measured after a 
pulse. Vc.o.m. will therefore offer a useful probe to study the 
impact of disorder on transport of strongly interacting atoms 
in optical lattices. 


III. DYNAMICS FROM EQUATIONS OF MOTION 


To derive the equations of motion we construct and solve 
the hierarchy of coupled differential equations with Hubbard’s 
decoupling. The commutator in Eq. (8) can be evaluated: 

= {lH>-m){ph) + U{n l ) 

+ Y [TijiPi'j) - (13) 

3 

where the U term contains a higher order correlator: 


To study the center-of-mass dynamics we derive equations 
of motion for correlation functions related to observables. The 
trapping potential in Eq. (1) breaks translational invariance. 
We will derive the equations of motion in the site (Wannier) 
basis as opposed to the more conventional fc-space (Bloch) ba¬ 
sis to allow studies of local dynamics in trapped lattices. We 
approximate the equations of motion by relying on the strong 
interaction/high temperature limit, Eqs. (5). The next section 
shows that our approximation reduces to Hubbard’s decou¬ 
pling of the equations of motion, the Hubbard-I approxima¬ 
tion [43], that introduced the concept of Hubbard-band quasi¬ 
particles. The equations of motion derived here therefore offer 
a tool to study the local dynamics of Hubbard-band quasipar¬ 
ticles in the absence of translational invariance. 

The exact equations of motion for the charge and spin de¬ 
grees of freedom are given by: 

ih Ml = {[ p a ii : H AH }) (8) 

and 

= ([Sr .z.^ah]), (9) 

respectively. Here the correlator: 


= Pinim-v-TH',-*)- (14) 


The central aim of our protocol is to numerically solve 
Eq. (13) and use the results to evaluate Eq. (6). This will re¬ 
quire an estimate for F^. 

To estimate F" ; we derive the equations of motion for 
this higher order correlation function as well. The oper¬ 
ator evolves as: dT^Jdr = dpfn/drfoi-tr — rn>_ a ) + 
pfq(dni t - a /dT — dnii -a / dr). We use this relation to ap¬ 
proximate the evolution of (T^) by inserting Eq. (13) and 
decoupling all products of Tf, L and pf n : 


^ = (pi- - fMi)( rfr) + U{V{, l ){ni- (T - ni’-o) 
+ (m-a - nv,-a) Y [ T l,j(Pl'j) ~ T l\j(Pjl)] 


+ ( Pin) Yj [ Ti p(Pi/ ~ Pji) ~ Ti'jiPi'j ~~ Pjv) 

3 


(15) 


The key decoupling used in deriving this equation is given by 
a Hartree-Fock-like decoupling in the equations of motion of 
the form: 


PljPvi -»• (Pif)Pl'i 

ni-aTfn -A (ni-a) (16) 


Pin = 4,* c i,* 0°) 

is the single-particle density matrix which is off-diagonal in 
the site indices, l and V, but measures density along the diag¬ 
onal, since p= np a . The spin operator is: Sjqj = ipterip,, 

where the fermion spinors are: ipj = (cjcj and cr are the 
Pauli matrices. The equations of motion can be generalized to 
include time dependence in the Hamiltonian but we exclude 
that case here. 

The high temperature limit studied here suppresses spin or¬ 
der (which emerges for kgT < t 2 /[/). This implies that for 
an equal number of atoms in each spin state we have a param- 
agnet: 

(pIi) = (pIi)- (ID 

To focus on the charge degrees of freedom deep in the para¬ 
magnetic regime we also assume an absence of in-plane spin 
order as well. This leads to: 

= (Sf,,) = (S?I i) = 0, (12) 

thus allowing us to focus on approximations to Eq. (8) only. 


The next section shows that this decoupling reduces to Hub¬ 
bard’s decoupling [43] that has been conventionally imple¬ 
mented in a Green’s function approach [43, 49], It is im¬ 
portant to note that this decoupling goes beyond conventional 
Hartree-Fock decouplings of the Hamiltonian [50] (which 
only capture the dynamics of very weakly interacting limits) 
to instead apply a decoupling in the equations of motion of 
higher order correlation functions. It therefore offers a ro¬ 
bust formalism that captures both weak (t U ) and strong 
(t <C U) interaction limits of the paramagnetic phase. 

We self-consistently solve Eqs. (13) and (15) for the time 
evolution of the correlation functions. We then use the cor¬ 
relation functions to evaluate the center-of-mass position and 
velocity. The time evolution of other local correlation func¬ 
tions can also be found. For example, the double occupancy, 
(ni^npi), can be obtained from: 

in<1 ^r = (w ' - w)<7?!) + {Pn)Y T uW) - <^T>] 

3 

+ ( ni-a) ^ [■ Tij(pi'j) — Ti r ,j(Pji)] 

3 

+ Wt)(i-(^w))(i-V) 


(17) 
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where the off-diagonal operator: 

7 VI = Phm,-*, (18) 

captures the conditional hopping of doublons. 


IV. CONNECTION TO HUBBARD’S APPROXIMATION 


In this section we argue that the formalism we have con¬ 
structed can be understood in a quasiparticle picture. In 
strongly interacting systems we often rely on mappings to 
weakly interacting quasiparticles to gain a quantitative un¬ 
derstanding of otherwise intractable problems. Quasiparticles 
therefore offer useful tools to probe many-body localization 
and related phenomena. We can then ask the following ques¬ 
tion that parallels inquiries into many-body localization of ele¬ 
mentary particles: Does spatial disorder localize weakly inter¬ 
acting quasiparticles at non-zero temperature? Here the inter¬ 
actions between the original particles are strong thus allow¬ 
ing significant dephasing from interactions. But quasiparti¬ 
cle problems are tractable and should therefore allow detailed 
quantitative studies. 

To connect the equations of motion to Hubbard-band quasi¬ 
particles we will show that our formalism reduces to Hub¬ 
bard’s approximation in the translationally invariant limit. 
Our formulation is a local theory designed to incorporate spa¬ 
tial inhomogeneity (i.e., trapping and disorder in the quasi¬ 
particle degrees of freedom). By assuming translational in¬ 
variance we can show that the above formalism simplifies to 
the equations of motion found from Hubbard’s approximation. 
We first briefly review Hubbard’s approximation. 

Hubbard’s approximation applies the Hartree-Fock decou¬ 
pling to the equations of motion for the Green’s functions. The 
approximation is, unlike the ordinary Hartree-Fock approxi¬ 
mation, exact in both the band limit (i.e., no interactions) and 
the atomic limit (i.e., infinite interactions). The approxima¬ 
tion assumes two Hubbard bands of quasiparticles where the 
band parameters are renormalized by the density and the in¬ 
teraction. Exact solutions of Hubbard’s equations of motion 
are possible in the translationally invariant limit (fij -A —//q 
in Eq. (2)). The quasiparticle Green’s function is: 


G k AE) 


h 

E-[e(k)-^ + E a (E)Y 


(19) 


where the nearest neighbor tunneling leads to the single parti¬ 
cle band dispersion: 


e(k) = —21 ^ cos (k u a). ( 20 ) 

is€x,y,z 


The self energy is [43]: 

£*(£) = 


U {n a )(E + p 0 ) 


( 21 ) 


E + p 0 - U(1 - (n-<j))' 

We can define a spectral density that is useful for calculations: 
S k , a {E) = HS[E - e{k) + p 0 - Z a (E )}. (22) 
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FIG. 2. (Color online) Plot of the energy versus wavevector for a 
translationally invariant lattice. The bottom (top) solid lines plot the 
energies of the lower (upper) Hubbard bands, -E(, >cr (fc), in one di¬ 
mension. The parameters U/t = 9 and (n) = 0.7 are chosen as 
characteristic of the center of system for the 7Er data in Table I 
(with u> = 0 and Vp = 0). The dotted line plots the non-interacting 
case, Eq. (20), but in one dimension, for comparison. Here we see 
that the lower Hubbard band is very similar to the non-interacting 
band. 


From the spectral density we find two (Hubbard) bands with 
spectral weights that depend on the density and interaction. 
The energies of each band are: 


Eb,*( k ) = 


U + e(k) 


+ (-l ) 6 


[u-s(k)Y 


U (n_ a )e(k), 


where b = 1 (b = 2) denotes the lower (upper) Hubbard band. 
In the limit of weak interaction the bands become degener¬ 
ate and we recover the Hartree-Fock limit from Hubbard’s ap¬ 
proximation. 

The Hubbard bands split in the strong interaction limit. To 
see this we expand Eb i<7 in powers of 1/U. We find: 

E ha = [1 -(n- a )}e(k) + 0(t 2 /U) 

E 2 ^ = U + (n_ a )£(k)+(D(t 2 /U) (23) 


This shows that, to lowest order, lower Hubbard-band quasi¬ 
particles can be thought of as non-interacting particles but 
with a renormalized hopping, t [1 — (n_ CT )]. (Technically, 
the renormalized hopping allows the quasiparticles to inter¬ 
act through the mean field) The upper Hubbard band is sim¬ 
ilar but with an energy offset U and a renormalized hopping 
t(n- a ). An important aspect of Eq. (23) is that the corrections 
are ~ t 2 /U and are therefore much smaller than the tempera¬ 
ture in most ongoing optical lattice experiments. Fig. 2 plots 
Eb,a{k) in comparison to s(k) along one dimension to show 
that the energetics of Hubbard-band quasiparticles in the low¬ 
est band are close to those of free particles. 

We now show that the formalism presented in Eqs. (13)- 
(17) reduces to Hubbard’s approximation in the translation- 
ally invariant limit. To show this we simplify the equations 
of motion for and 7 ^. We can then solve the equations 
of motion by Fourier transforming into energy and wavevec¬ 
tor variables. We find that the resulting energies are given by 
Eb,a{k). 
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Eqs. (13) and (17) define a coupled set of equations that can 
be solved analytically in the translationally invariant limit. We 
note that these equations are coupled since = 7 — 'yfj. 
We impose translational invariance by setting pi = /i/<. The 
density then becomes uniform: (n CT ) = {np a ). We Fourier 
transform all terms in Eqs. (13) and (17). For example, we 
set: 

Ply = N- 1 J2 (24) 

i,v 

where N s is the number of sites. 

We can also transform the coupled set of first order differ¬ 
ential equations in time to a set of coupled algebraic equations 
by Fourier transforming to energy space. We then find: 

- Epl k , = U{r k y 7# J + H-fc) - <k')\ Ply 

-Elly = Elly + (n-a) [e(-fe) - e(fc')] Ply> ( 25 ) 

where we have dropped higher order terms, i.e., terms of 
the form U{n- a ) ll k f. We have also made use of 2];/ = 

N7 1 E fe e(fc)e lfc ' (Ri ’ _Ri,) . 

Eqs. (25) can be solved analytically for the eigenvalues E 
by setting k! = 0 and solving for p k 0 and 7 ^ 9 . We can, 
without loss of generality, set e(0) = 0 in Eq. (25) to make 
contact with the Hubbard approximation. We find three dis¬ 
tinct modes. A trivial high energy mode with E = U corre¬ 
sponds to a non-dispersive doublon mode obtained from so¬ 
lutions with p k0 = 0. But the remaining two modes we find 
have precisely the same energies as those found in Hubbard’s 
approximation: £), jCr (fe). We have therefore shown that the 
formalism presented in Eqs. (13)-(17) reduces to Hubbard’s 
approximation in the translationally invariant limit. 

The reduction of the transport problem posed by Eq. (1) in 
a high temperature paramagnetic limit into that of transport of 
Hubbard-band quasiparticles has two important implications. 
The first is practical: We will be able to compute correlation 
functions for the initial state using the spectral density. This 
is discussed in Sec. V. 

The second implication is phenomenological. Since the 
strongly interacting limit can be thought of as nearly free 
quasiparticles, the addition of disorder should show features 
qualitatively similar to a weakly interacting system. We have 
verified that the quasiparticle picture remains valid even for 
large disorder strengths, As/i ~ 40, by checking that the 
Hubbard-band spectral weight is non-zero. We can therefore 
construct an effective model of quasiparticles in a disordered 
lattice (but in the absence of a trap): 

Hefl = beT + Pj,b,anj,b,a (26) 

j,j',cr,b j,b,<7 

where the tilde indicates quasiparticle operators, p, is the 
chemical potential renormalized by the self energy and Tjji 
indicates quasiparticle hopping with: 

T)'° = N- 1 ]T S b , CT (fc)e ifc (Rl - R '' ) . (27) 

fc 



FIG. 3. (Color online) Schematic of one two-dimensional plane of 
the site-disordered cubic lattice. The blue (red) spheres depict spin 
up (down) quasiparticles in the lowest Hubbard band. The high¬ 
lighted nearest neighbor bonds between particles of opposite spin 
symbolize the renormalized quasiparticle hopping in Eq. (26). The 
hopping is suppressed on average. For example, a lower-band quasi¬ 
particle with spin a has a renormalize nearest neighbor hopping 
t(l — n —er), for large U. 


Here we have assumed that the quasiparticle energies, 
Eb,a{k), depend on the Fourier transform of the randomly 
distributed density. 

We can get an intuition for the renormalized hopping if we 
assume that the density, on average, remains uniform in the 
presence of disorder. Eqs. (23) show that in the strongly in¬ 
teracting limit this renormalized hopping reduces to: T- ^ « 

T jtj/ (l-n^)+0(t 2 /U) and iff, « T jtj > (n_ CT ) +0(t^/U), 

for the lower and upper Hubbard bands, respectively. The 
renormalized hopping is shown schematically in Fig. 3. 

H e ff is an effective theory of quasiparticles that must, in 
principle, be solved self-consistently. But it should nonethe¬ 
less reveal a mobility edge of Hubbard-band quasiparticles 
because it is essentially a non-interacting Anderson model of 
Hubbard-band quasiparticles. For example, it is well known 
that the Anderson model in the cubic lattice with uniform 
disorder exhibits a mobility edge at A u/B = X c , where 
X c « 1.6 [51], For the lower Hubbard band in the absence of 
a trap and in a paramagnetic state the quasiparticle bandwidth 
becomes B = 12f(l — (n)/ 2). iT e ff therefore qualitatively 
predicts a mobility edge for Hubbard-band quasiparticles. We 
will return to this point in discussing the suppression of trans¬ 
port in Sec. VIII. 

iT e ff a ls° shows that ignoring disorder in / and U is justified 
in the large U limit. It Ref. [30] it was shown that speckle 
disorder leads to a narrow Lorentzian distribution of / and U. 
Even though the distribution is narrow, these parameters could 
in principle have significant contributions to transport due to 
the tails of the distribution. But we note that the large U limit 
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is dominated by transport of quasiparticles (not the original 
particles). Eqs. (23) and (26) explicitly show that the effec¬ 
tive quasiparticle hopping, T, and chemical potential, p, are 
implicitly disordered even if disorder in t and U are excluded. 
This shows that excluding disorder in t and U still leaves an 
effective model with all terms disordered. Including disorder 
in t and U should therefore not qualitatively impact the trans¬ 
port properties of quasiparticles in the large U limit. 


Using the these relations we are able to set the initial state 
correlation functions with a protocol discussed in Sec. VII. 
The protocol allows use of the Hubbard approximation to 
compute initial state correlation functions at fixed entropy for 
a given disorder configuration. The following discusses the 
temperature dependence in the initial state in the presence of 
disorder. 


V. INITIAL STATE 


VI. ADIABATIC HEATING DUE TO DISORDER IN THE 
INITIAL STATE 


To time evolve correlation functions we must accurately es¬ 
tablish the initial state. The system evolves in the absence of 
a heat or particle number bath. The short-time dependence 
therefore crucially depends on the initial state. We note that 
the Hubbard approximation is very accurate in the limit de¬ 
fined by Eqs. (5). To see this we note that the static properties 
of optical lattice experiments with two-component fermions 
are also accurately captured by a high temperature series ex¬ 
pansion of Eq. (1) [52-54], We have checked that the high 
temperature series expansion and the Hubbard approximation 
agree in the limits discussed here, Eqs. (5). 

In this high T regime the initial state is also accurately cap¬ 
tured by the local density approximation [49, 53], We take 
each site as a uniform system and compute correlation func¬ 
tions. In the local density approximation we assume that each 
site of the trapped system can be approximated with parame¬ 
ters for a uniform system by setting /// to be the chemical po¬ 
tential for the I th uniform system, and we average over all uni¬ 
form systems. Correlation functions from each site are then 
combined. In the case of multi-site correlation functions a 
complication arises: the chemical potential varies from site to 
site. Here we find that nearest neighbor correlation functions 
are sufficient to describe the initial state, since long range cor¬ 
relation functions decay quickly at these temperatures. As a 
result, we are able to approximate two-site correlation func¬ 
tions by setting the chemical potential to be the average be¬ 
tween the neighbors. 

The initial state can be approximated using correlation 
functions computed directly from the spectral function within 
the local density approximation. The spectral theorem implies 
that we can compute the initial (r = 0) correlation function 
for pin using: 


(Pi'i)(t = °) 


g—i(H;/ —Ri)k 

^ 2 HN S 

k 

Sk,a(E - p), 


dEf(E) 


(28) 


where p = (/q — /q/ )/2 and f(E) is the Fermi-Dirac distribu¬ 
tion function. Here we assume (pfn) is equal to its Hermitian 
conjugate. A similar relation can be used to obtain as 
well: 


{Pi'i n l,-a)(r = 0 ) = 


0 -i(R lt -Rl)-k r°° 


2HN, 


dEf(E) 


E — e(k) 
U 


S k , a {E~p). (29) 


The temperature in ultracold atom experiments is deter¬ 
mined by the entropy. The relationship between temperature 
and entropy relies, in general, on the intricate interplay be¬ 
tween kinetics and interactions. The addition of disorder adds 
another complication that alters the entropy-temperature re¬ 
lation. Below we show that the addition of disorder leads 
to adiabatic heating in the initial state. Specifically we find 
that, at fixed entropy, increasing disorder increases the tem¬ 
perature. This observation has important consequences for 
the interpretation of the data in Ref. [39] and other optical 
lattice experiments because increasing disorder strengths also 
increases temperature. In subsequent sections we take adia¬ 
batic heating from disorder into account when preparing the 
initial state in a trap. 

We use the high temperature series expansion to show that 
the paramagnet experiences adiabatic heating due to disorder. 
The solid line in Fig. 4 shows an example of the entropy per 
particle versus temperature for an initial state without a trap. 
We set po/t = 3.8 because it characterizes the non-disordered 
limit of experiments reported in Ref [39]. We find (n) < 1. 
Here we see that a fixed entropy (horizontal dotted line) sets 
a low temperature, Tq, in the absence of disorder. Because 
optical lattice experiments take place in the absence of a heat 
bath, entropy is preserved when a disordered optical lattice is 
applied to a trapped gas. We then include a disorder strength 
A(/ = 20f in a calculation of the entropy per particle. We use 
the local density approximation and integrate over disorder 
configurations (See Eqs. (A3)). The dashed line shows the 
disorder averaged results. The entropy is significantly lower. 
The system therefore acquires a higher temperature, 7)/, at 
the same entropy. 

Adiabatic heating due to disorder arises because increas¬ 
ing the disorder strength in a single band reduces the num¬ 
ber of available states. As a result the entropy (which is the 
logarithm of the number of available states) decreases with 
increasing disorder. The net effect is then an increase in tem¬ 
perature if the entropy is required to be fixed while increasing 
disorder. 

Adiabatic heating becomes more pronounced with expo¬ 
nentially distributed disorder. Fig. 5 plots the entropy per par¬ 
ticle as a function of both exponential disorder strength and 
temperature. The black lines depict adiabats. The correspond¬ 
ing temperature can therefore increase by as much as a factor 
of 2 at fixed entropy over the range of disorder strengths con¬ 
sidered here. The impact of adiabatic heating due to disorder 
on center-of-mass transport in trapped systems is discussed in 
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FIG. 4. (Color online) The disorder-averaged entropy per particle 
computed as a function of temperature for Eq. ( 1 ) in the absence of 
a trapping potential or a pulse (to = 0 and Vp = 0). The horizontal 
dotted line indicates a fixed entropy per particle, S/N = 1.2fcB- The 
solid (dashed) lines were computed using Ap = 0 (Ay = 20f) and 
fMi/t = 3.8. The vertical lines labeled with Ty and Th point to low 
and high temperatures, respectively. The entropy-temperature curve 
with a high disorder leads to a higher temperature. 


more detail in the following sections. 


VII. CENTER-OF-MASS DYNAMICS: COMPARISON 
WITH EXPERIMENT 

This section culminates in a direct comparison between re¬ 
sults from the equations of motion and experiments. We find 
that small system size simulations can be scaled to directly 
compare with experiments with no fitting parameters. The 
close comparison between experiment and theory shows that 
we can interpret the experiments of Ref. [39] as transport of 
Hubbard-band quasiparticles. The simulations and experi¬ 
ments are consistent with finite size precursors of Anderson 
localization of Hubbard-band quasiparticles. 

We now use Eqs. (28) and (29) to compare with experi¬ 
ments in Ref. [39] using experimental input parameters from 
Table I. To use our formalism to compute the center-of-mass 
dynamics we prepare an initial state at fixed entropy in a dis- 



FIG. 5. (Color online) Disorder averaged entropy per particle plotted 
as function of both temperature and disorder strength for Eq. (1) in 
the absence of a trapping potential or an applied force (i.e., to = 0 
and Vp = 0). The 8 th order high temperature series expansion was 
used within the local density approximation. U/t = 9 and p.o/t = 
3.8 were chosen as characteristic of the center of system for the 7Ep 
data in Table I. The black contour lines indicate adiabats that reveal 
significant adiabatic heating due to increasing exponential disorder. 


ordered landscape. The system is numerically time evolved. 
The center-of-mass velocity is computed at the pulse time and 
then disorder averaged. These simulations are performed on 
system sizes up to L = 11, with L = L x = L y = L z . 
Finite size extrapolation is performed by decreasing the trap 
frequency and repeating the simulation for large system sizes 
while keeping // 0 fixed to values found for experimentally rel¬ 
evant system parameters. 

To keep the pulse time short on the time scales of the trap¬ 
ping frequency (as is done experimentally [39]) we have to 
rescale the pulse time used in our simulations. The pulse time 
at system size L, Ty, is adjusted for each trap frequency at 
system size L, wy, to maintain Ty = Tp^/w/wy. This al¬ 
lows a scaling to the trapping frequency and the pulse time 
found in Table I, u> and rp, respectively. The impulse formula 
(See Sec. A) shows that this establishes an uj l scaling of 
Vc.o.m.- This scaling is expected since the center-of-mass ve- 

_ 1 /2 

locity from the impulse formula scales as Vo ~ ry ~ 

(see Sec. A). We have checked below that our finite size ex¬ 
trapolations do scale as to L 1 /2 , as expected. 

We use the following protocol to prepare initial states: 1) 
We choose an entropy per particle determined experimentally, 
high trap frequency (chosen to trap the system within the fi¬ 
nite size limitations of our simulations), and a small number 
of particles. 2) We choose a random distribution of chemical 
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FIG. 6. (Color online) Top: Disorder averaged center-of-mass veloc¬ 
ity as a function of the disorder strength for two different entropies. 
Here the initial state correlation functions are estimated in the lo¬ 
cal density approximation (using Eqs. (28) and (29) in combination 
with a high temperature series expansion) and time evolved in the 
trap (using Eqs. (13) and (15)). Model parameters are taken from the 
6 Er data in Table I. The S/N = 1.5kB results are plotted only for 
large disorder strengths because here adiabatic heating allows access 
to temperatures high enough to be consistent with the approxima¬ 
tions made in preparing the initial state. Bottom: The circles plot 
the same as the top panel and the diamonds plot experimental data 
from Ref. [39] for comparison. The lines are a guide to the eye. The 
error bars on the numerical simulations are the standard error found 
from disorder averaging, while the experimental error bars are the 
standard error in the mean for 7-9 measurements averaged for each 
point. 


potentials according to Eq. (3). 3) We then self-consistently 
adjust fio and T so that the particle number and entropy match 
the values set in step 1. This is done using a high temperature 
series expansion in the local density approximation. The se¬ 
ries expansion is controlled a these temperatures because we 
can check higher orders [53, 54]. We find that 8 th order in the 
expansion is sufficient for parameters considered here. The 
Hubbard approximation gives identical results for thermody¬ 
namic functions. 4) We then use Eqs. (28) and (29) to compute 
the initial state correlation functions. 5) We then return to step 
1 to repeat the process with a smaller trap frequency. 

We find that adiabatic heating in the initial state increases 
the temperature by no more than a factor of 2. For all sys¬ 
tem sizes studied we find that the temperature remains nearly 
constant as function of system size. At the largest disorder 


FIG. 7. (Color online) The same as Fig. 6 but for the 7Er data in 
Table I. Here the comparison between theory and experiment is better 
because U/t is larger. 


strengths, A# ~ 1.5 Er, we still find UrT < At. We con¬ 
clude that adiabatic heating increases the temperature but the 
temperature is still well below the bandwidth, 12i. 

Given the initial state, we numerically time evolve correla¬ 
tion functions according to Eqs. (13) and Eqs. (15), extrapo¬ 
late to the thermodynamic limit, and disorder average. Figs. 6 
and 7 plot Vc.o.m. versus disorder strength for the 6 Er and 
the 7 Er parameters, respectively. The data result from time 
evolving the initial correlators, Eqs. (28) and (29). The top 
panels show results for two different entropies. The larger en¬ 
tropy leads to temperatures with T > t. The approximations 
made here (paramagnetic order, no spin correlations, and the 
local density approximation) are therefore valid at all disorder 
strengths for the higher entropy. The top panels also compare 
low entropy data that is consistent with the entropies used in 
experiments (see Table I). Here adiabatic heating increases the 
temperature to T > t only for A r > 0.2 Er. Below these 
disorder strengths the approximations made here break down 
because the temperatures are low enough to introduce poles 
in thermodynamic functions using either the high temperature 
series expansion (even out to 10 th order) or the Hubbard ap¬ 
proximation. 

The top panels of Figs. 6 and 7 clearly show a suppres¬ 
sion of the center-of-mass velocity with disorder. The map¬ 
ping to quasiparticles in the lowest Hubbard band allows de¬ 
lineation of the sources of the suppression: 1) As exponen¬ 
tially distributed disorder is increased, the bias in the distri- 
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bution leads to more sites with higher densities. The increase 
in average density slows the propagation of the Hubbard-band 
quasiparticles because the renormalized tunneling is given by 
i(l — n/2). This effect was implicit in the suppression shown 
in Sec. A (See Fig. 8). We find that this is a weak effect be¬ 
cause the system is dilute, i.e., (rij) / 2 <C 1 for many sites (the 
edges make up about 1/3 of the system) 2) Adiabatic heat¬ 
ing due to disorder also suppresses Vc.o.m.- The increase in 
the resulting temperature lowers the nearest neighbor corre¬ 
lations, e.g., (pfi + i), inherent in the initial state. The ini¬ 
tial state is therefore slower to respond because Vq .o.m. scales 
linearly with terms like t(pf l+1 ). This effect was shown to 
dominate only at lower disorder in Sec. B (See Fig. 9). Fur¬ 
thermore, we find that the temperature is at most B /3 at the 
largest disorder strength. A# ~ 1.5-E/j. 3) These effects are 
modest and are not sufficient to completely localize the center 
of mass. The final effect derives from disorder induced scat¬ 
tering. The presence of disorder lowers the localization length 
so that propagation is impossible for Ae > 0.5 Eft,. This fi¬ 
nal effect is consistent with a finite size precursor of Anderson 
localization of Hubbard-band quasiparticles because the crit¬ 
ical disorder strength, A e ~ 0.5 Er, is near the approximate 
location expected for the Anderson metal-insulator transition, 
near B « 0A7Er. 

The bottom panels in Figs. 6 and 7 show a comparison be¬ 
tween the results obtained from our formalism and the ex¬ 
perimental data of Ref. [39]. The comparison is made where 
possible (in the high temperature regime). The agreement in 
Fig. 7 is better because U is larger. The Hubbard approxi¬ 
mation is technically a strong coupling approximation that is 
exact in the U/t —> oo limit. The comparison suggests that the 
data from Ref. [39] can be thought of as revealing a mobility 
edge of Hubbard-band quasiparticles. 


VIII. DISCUSSION 

We have found that two-component fermions in an optical 
lattice fail to respond to a force and undergo mass transport 
for sufficiently strong disorder, implying a phenomenon remi¬ 
niscent of Anderson localization in bulk systems. At strong 
disorder strengths the atoms fail to move under weak per¬ 
turbations. Here the suppression of quantum diffusion indi¬ 
cates that the assumption of a thermal initial state is incor¬ 
rect, i.e., that the system is inherently non-ergodic at large 
disorder strengths. Our comparison between theory and ex¬ 
periment are therefore consistent with Anderson localization 
of Hubbard-band quasiparticles at large disorder strengths but 
a mobile state of Hubbard-band quasiparticles at low disorder 
strengths. We interpret these results as evidence for a mobility 
edge of Hubbard-band quasiparticles. 

We can compare the center-of-mass velocity studied here 
with conductivity studied in solids. Both measures can be 
used as diagnostics of localization. The DC Conductivity 
in solids is typically defined in infinite system sizes. The 
DC conductivity therefore gives a long time/large length scale 
probe of the single particle density matrix. The center-of-mass 
velocity is proportional to mobility and therefore also offers 


an equivalent probe of the single particle density matrix pro¬ 
vided the system is infinitely large and it is allowed to evolve 
indefinitely. But the center-of-mass velocity studied here was 
considered on time scales inversely proportional to the trap 
frequency and in finite system sizes. We therefore conclude 
that the results presented in Figs. 6 and 7 only offer a finite 
size estimate for the conductivity. 

Our work opens interesting directions for future studies of 
localization physics with Hubbard-band quasiparticles. The 
work presented here is consistent with quantum Monte Carlo 
results [16] and dynamical mean field theory studies of the 
Anderson-Hubbard model [14, 17, 18]. But these methods 
could be used to tackle lower temperature limits and include 
spin fluctuations in a comparison with low temperature exper¬ 
iments. 

Furthermore, future work will be needed to rigorously es¬ 
tablish a connection between the localized state found here 
and many-body localization. The suppression of transport at 
non-zero temperatures found here is a necessary condition for 
many-body localization. But future work should look at suffi¬ 
cient conditions for many-body localization using, e.g., entan¬ 
glement measures in the Anderson-Hubbard model, to make a 
direct comparison with experiments. 

In preparing this manuscript we became aware of work in 
Ref. [55] that compared the entanglement entropy with popu¬ 
lation imbalance in incommensurate optical lattices. 

V.W.S. acknowledges support from AFOSR under grant 
FA9550-11-1-0313. B.D. acknowledges support from the 
NSF under grant PHY12-05548 and from the ARO under 
grant W9112-1-0462. 

Appendix A: Order of Magnitude Estimate 

This section uses a semiclassical impulse formula for 
Hubbard-band quasiparticles to estimate the center-of-mass 
velocity dependence on disorder strength for very weak dis¬ 
order. This estimate shows that renormalization of the quasi¬ 
particle hopping due to disorder can suppress the center-of- 
mass velocity. It also yields the correct order of magnitude 
for the center-of-mass velocity at low disorder. A simple or¬ 
der of magnitude estimate for the center-of-mass velocity will 
be useful in establishing a scaling relation to extrapolate our 
finite sized simulations to experimental system sizes. 

To estimate the center-of-mass velocity we use a semi¬ 
classical estimate of velocities in combination with the local 
density and effective mass approximations. The quasiparti¬ 
cle effective mass in the lowest Hubbard band is obtained 
from the single-particle effective mass using the replacement 
t —► t(l — rij/2 ): 

m *> = 2i(l — ((nj)) D /2)a 2 ’ (A1) 

where the limit (rij) -A 0 returns the single-particle effective 
mass. Note that disorder averaging is implicit in this defini¬ 
tion. 

At short times, the semiclassical estimate of the center-of- 
mass velocity reduces to the well known impulse formula. We 



11 


apply the impulse formula to the dynamics of Hubbard-band 
quasiparticles in the lowest band. (Note that the impulse for¬ 
mula also follows from the generalized Kohn’s theorem in 
an effective mass approximation) Averaging the velocity of 
each site, (f?J), leads to a total center-of-mass velocity for 

one disorder configuration TV -1 ZZf s (nj)(R *). Applying the 
impulse formula to Hubbard-band quasiparticles and averag¬ 
ing over disorder realizations gives an approximation to the 
center-of-mass velocity: 

N s 

v 1 = Vo[l - ^((n J -) 2 )c/(2Af)], (A2) 

j 

where Vo = 2aVpTpt/H 2 depends linearly on rp and 
(( n j) 2 )D indicates the disorder average of (nj) 2 . 

V 1 gives the correct order of magnitude for the center-of- 
mass velocity. To show this we use the high temperature series 
expansion to estimate the density in the initial state in the trap. 
We choose the parameters for the 7 Ep lattice depth presented 
in Table I but we fix the entropy to be S = 1.9 kp- 

We use a simplified version of the protocol constructed in 
the main text to get a rough estimate of V 1 . Once the entropy 
and particle number are fixed, the approach used in the main 
text then finds the //(, and T at each disorder configuration 
using the high temperature series expansion. These parame¬ 
ters are then, for each disorder configuration, used to compute 
( rij ) within the trap. Disorder averaging proceeds by sum¬ 
ming the center-of-mass velocity over all disorder configura¬ 
tions. But in this section we solve for the chemical potential 
and temperature differently so we can access experimentally 
relevant system sizes without finite size extrapolation. We use 
the high temperature series expansion to approximate the en¬ 
tropy and density with integration (rather than explicit sum¬ 
mation) over the disorder distribution: 

nOO 

((S))d~ / deP E (e)S(e) 

Jo 

/»oo 

(( n)) D ~ / deP E (e)(n(e)). (A3) 

Jo 

These approximations can be used to self-consistently solve 
for T and // (1 given S and N for large systems sizes. This 
simplified protocol uses entropies and densities that are not 
self-consistently solved for each disorder configuration but are 
instead taken in a mean-field limit separately. As a result, 
self-consistent solutions of these coupled formulas only offer 
a rough estimate for T and /i () because they are assumed to 
decouple for each disorder configuration. We can therefore 
only apply these approximations for low disorder strengths. 

The top panel of Fig. 8 plots the disorder-averaged density 
of the central site in the trap as a function of disorder strength. 
Here we see that the density decreases due to adiabatic heat¬ 
ing and a redistribution of the particles due to biased exponen¬ 
tial disorder. The quasiparticle effective mass (middle panel) 
therefore also decreases. 

The bottom panel of Fig. 8 plots the disorder-averaged 
center-of-mass velocity from Eq. (A2). Here we see that the 
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FIG. 8. (Color online) The top and middle panels plot the density 
and the quasiparticle effective mass, Eq. (A1 ), respectively, as a func¬ 
tion of disorder strength for a chemical potential at the central site. 
The bottom panel plots the center-of-mass velocity versus disorder 
strength from a disorder-averaged impulse formula, Eq. (A2), that es¬ 
timates the velocity of Hubbard-band quasiparticles in the trap size 
consistent with experiment. The local density approximation was 
used to sum over all sites. All quantities are computed using the high 
temperature series expansion at 8 th order with the parameters cho¬ 
sen from the 7Er data in Table I but with S/N = 1.9 k B . Eqs. (A3) 
were used as rough estimates for disorder averaging. 


velocity decreases due to an enhancement of the density. The 
experimental data, for comparison, starts out with a center-of- 
mass velocity ~ 1 mm/s. The impulse formula for Hubbard- 
band quasiparticles therefore gives the correct order of magni¬ 
tude and shows suppression due to a modulation of the density 
due to disorder. 


Appendix B: Temperature Dependence 

In this section we study the temperature dependence of the 
center-of-mass velocity in small trapped systems by solving 
for the dynamics of correlators using Eqs. (13) and (15). Here 
it is shown that temperature increases (expected in adiabatic 
heating) suppress the center-of-mass velocity but only for low 
disorder strengths. 
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FIG. 9. (Color online) The disorder averaged center-of-mass velocity 
(Eq. (7)) as a function of temperature for several disorder strengths 
computed from solutions of Eqs. (13) and (15). Parameters are cho¬ 
sen to yield a small system size replica of the larger system implied 
by the parameters for the 7 Er data in Table I (see text). The veloc¬ 
ities are disorder averaged using uniform disorder. Here we see that 
increasing temperature suppresses the velocity only at low disorder 
strengths. 


We can use Eqs. (13) and (15) to compute the center-of- 


mass dynamics in trapped systems on small system sizes. We 
solve Eqs. (13) and (15) numerically. The initial state is de¬ 
termined using Eqs. (28) and (29) within the local density ap¬ 
proximation at fixed temperature. Fig. 9 shows example re¬ 
sults for the center-of-mass velocity. The simulations are car¬ 
ried out on a periodic cubic lattice with edges of size L = 11 
where the trap zeroes the density at the the edges. We consider 
a small system size replica of larger experimental parameters 
by choosing a stronger trap frequency hu>/t = 0.757 but at 
the same chemical potential as that found for experimental 
system sizes, fi 0 /t = 3.8. rp = 0.514 h/t is chosen by a 
trap-dependent rescaling discussed in Sec. VII. The entropy 
is allowed to vary but otherwise the remaining parameters are 
chosen from the 7 Er data in Table I. 

Fig. 9 shows that by increasing temperature, the center-of- 
mass velocity can decrease at low disorder. This is the oppo¬ 
site of what is expected from variable range hopping in com¬ 
mon regimes, e.g., in semiconductors, where the presence of 
a bath typically increases conductivity with increasing tem¬ 
perature. Here we do not have an external bath. At low 
disorder, increasing temperature suppresses the amplitude for 
particles to tunnel between neighboring sites, e.g., t(pf l+1 ), 
in the initial state. As a result, the center-of-mass velocity 
(which scales linearly with the nearest neighbor elements of 
the single-particle density matrix) is suppressed with increas¬ 
ing temperature. The high disorder limit has a different be¬ 
havior. Here the dynamics is strongly suppressed by disor¬ 
der and the thermal suppression of tunneling has little effect. 
These qualitative trends show that, when we study the experi¬ 
mentally relevant fixed entropy case, adiabatic heating due to 
disorder will tend to suppress the center-of-mass velocity only 
at low disorder strengths. 
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